

%% load theta from last iteration
load('parameter_estimates_direction_round1_replication.mat')
[~,use] = min(l0);
thetaStar = theta(:,use); %#ok<*NASGU>
clearvars -except thetaStar specif


%% load everything of interest
load('sim_eq_direction_round2.mat')
load('data_direction.mat')
[DMES,intervEq] = createDMES(intervEq,V_R,V_US,V_UK,V_France,V_Russia,V_China,YY,specif);


%% normalize continuous variables
data.gdp_pc = (data.gdp_pc-mean(data.gdp_pc))/std(data.gdp_pc);
data.pop = (data.pop-mean(data.pop))/std(data.pop);
data.polity = (data.polity-mean(data.polity))/std(data.polity);
data.terrain = (data.terrain-mean(data.terrain))/std(data.terrain);
dd = [data.USdist,data.UKdist,data.FRdist,data.RUdist,data.CHdist];
dd = (dd-mean(dd(:)))/std(dd(:));
data.USdist = dd(:,1);
data.UKdist = dd(:,2);
data.FRdist = dd(:,3);
data.RUdist = dd(:,4);
data.CHdist = dd(:,5);
clear dd


%% estimates of mean utilities
X=[ones(N,1) data.terrain data.gdp_pc data.polity data.pop]; K=size(X,2);
Z_US=[data.USdist data.USally data.UScolony data.USwar]; KZ=size(Z_US,2);
Z_UK=[data.UKdist data.UKally data.UKcolony data.UKwar];
Z_France=[data.FRdist data.FRally data.FRcolony data.FRwar];
Z_Russia=[data.RUdist data.RUally data.RUcolony data.RUwar];
Z_China=[data.CHdist data.CHally data.CHcolony data.CHwar];
Z=[Z_US,Z_UK,Z_France,Z_Russia,Z_China];
reb={1:K,2:KZ,2:KZ};
mpow={[1,3:K],[1,3:K],1:KZ,1:KZ};
Y=arrayfun(@(n)find(sum(YY==table2array(outcomes(n,4:end)),2)==size(YY,2)),1:N)';

options=optimset('Display','off','HessFcn',...
    @(x,lambda)HessPhat(x,lambda,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow));

% "long" replication
% thetaStarExp = [thetaStar(1:2);zeros(3,1);thetaStar(3:7);zeros(3,1);thetaStar(8:12);zeros(3,1);...
%     thetaStar(13:20);thetaStar(21:28);thetaStar(29);zeros(3,1);...
%     thetaStar(30);zeros(3,1);thetaStar(31:end)]; % expand thetaStar
% load('parameter_estimates_direction_round2_replication.mat','SV','srngSV')
% rng(srngSV)
% theta0=thetaStarExp+randn(length(thetaStarExp),SV)/10; % starting values
% theta=theta0; 
% l0=1:SV; % maximized likelihood values across starting values
% exitn=1:SV; % Knitro exit flags across starting values
% f0 = exitn;
% time=0;
% parfor sv=1:SV
%     tic
%     f0(sv) = Phat(theta0(:,sv),Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow);
%     if ~isfinite(f0(sv))
%         theta(:,sv) = NaN(length(theta0(:,sv)),1);
%         exitn(sv) = NaN;
%         l0(sv) = NaN; 
%         time=time+toc;
%     else
%         [theta(:,sv),l0(sv),exitn(sv)]=knitromatlab(@(x)Phat(x,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow),...
%             theta0(:,sv),[],[],[],[],[],[],[],[],options,'knitro1.opt'); % use knitro0.opt to check derivatives
%         time=time+toc;
%     end
% end

% "short" replication
load('parameter_estimates_direction_round2_replication.mat')
[~,use] = min(l0);
time=0;
[theta(:,use),l0(use),exitn(use)]=knitromatlab(@(x)Phat(x,Y,X,Z,YY,V_R,V_US,V_UK,V_France,V_Russia,V_China,U0,intervEq,DMES,reb,mpow),...
    theta(:,use),[],[],[],[],[],[],[],[],options,'knitro1.opt'); % use knitro0.opt to check derivatives
time=time+toc;


%%
clearvars -except exitn K l0 N SV theta theta0 time X Y YY Z reb mpow srngSV specif
save('parameter_estimates_direction_round2.mat')


%% verify optimality
[~,use] = min(l0);
thetaStar = theta(:,use);
load('parameter_estimates_direction_round2_replication.mat')
[~,use] = min(l0);
thetaStar0 = theta(:,use);
disp(' ')
disp('verify round 2 direction results:')
disp(['max. rel. diff. between estimates = ',num2str(max(abs(thetaStar0-thetaStar)./abs(thetaStar0)))])



